Read Bismark files
Negative control
## Reading ./input/NC_rep1.cov.gz
## Reading ./input/NC_rep2.cov.gz
## Reading ./input/NC_rep3.cov.gz
## Reading ./input/NC_rep4.cov.gz
## Reading ./input/NC_rep5.cov.gz
## Reading ./input/NC_rep6.cov.gz
## Hashing ...
## Collating counts ...
Simulated data
## Reading ./input/sim_rep1.cov.gz
## Reading ./input/sim_rep2.cov.gz
## Reading ./input/sim_rep3.cov.gz
## Reading ./input/sim_rep4.cov.gz
## Reading ./input/sim_rep5.cov.gz
## Reading ./input/sim_rep6.cov.gz
## Hashing ...
## Collating counts ...
Filtering and plotting
Function
beta_calc <- function(objectDGE, coverageThreshold, anno = anno) {
## Function to calculate Beta values from DGE object obtained from readBismark2DGE function.
## objectDGE: DGEList
## coverageThreshold: Coverage threshold to be used for all samples
## Return: Coverage, Coverage melted table, distribution, Beta value table, Beta value melted table, distribution
# Counts (raw coverage)
df <- data.frame(objectDGE$counts, stringsAsFactors = F, check.names = F)
me.raw <- df[, grep("Me", colnames(df))]
un.raw <- df[, grep("Un", colnames(df))]
cov.raw <- me.raw + un.raw
# Filter based on coverage
data <- NULL
if (is.null(coverageThreshold)) {
data <- df
coverageThreshold <- "NO filtering"
} else {
data <- df[rowSums(cov.raw > coverageThreshold) == ncol(cov.raw), ]
}
# Beta valus after filtering
me <- data[, grep("Me", colnames(data))]
un <- data[, grep("Un", colnames(data))]
beta <- me / (me + un)
colnames(beta) <- as.character(sapply(colnames(beta), function(x) strsplit(x, "-")[[1]][1]))
# Beta values for plotting
reshape.beta <- data.frame(t(beta))
reshape.beta$Samples <- rownames(reshape.beta)
reshape.beta.melt <- melt(reshape.beta)
# Plotting Beta value
p1 <- ggplotly(ggplot(reshape.beta.melt, aes(x = value)) +
geom_histogram(aes(y = ..density..), alpha = 0.7, fill = "#333333") +
geom_density(fill = "#CC6666", alpha = 0.3) +
theme(panel.background = element_rect(fill = "#ffffff")) +
facet_wrap(~Samples, nrow = 3, shrink = T) +
ggtitle(paste0("Beta distribution (coverage > ", coverageThreshold, ")")) +
labs(x = "Beta value", y = ""))
# Coverage plot after filtering
cov <- me + un
colnames(cov) <- as.character(sapply(colnames(cov), function(x) strsplit(x, "-")[[1]][1]))
reshape.cov <- data.frame(t(cov))
reshape.cov$Samples <- rownames(reshape.cov)
reshape.cov.melt <- melt(reshape.cov)
reshape.cov.melt.log2 <- reshape.cov.melt
reshape.cov.melt.log2$value <- log2(reshape.cov.melt.log2$value)
p2 <- ggplotly(ggplot(reshape.cov.melt.log2, aes(x = value)) +
geom_histogram(aes(y = ..density..), alpha = 0.7, fill = "#333333") +
geom_density(fill = "#20AAEA", alpha = 0.3) +
theme(panel.background = element_rect(fill = "#ffffff")) +
facet_wrap(~Samples, nrow = 3, shrink = T) +
ggtitle(paste0("Coverage distribution (coverage > ", coverageThreshold, ")")) +
labs(x = "Coverage (log2 scaled)", y = ""))
# PCA plot
## Coverage
pca.cov <- prcomp(t(cov))
p3.cov <- plot_ly(
x = pca.cov$x[, 1], y = pca.cov$x[, 2],
data = anno, color = ~Group, mode = "markers",
symbol = ~ as.numeric(as.factor(Group)),
marker = list(size = 10), hoverinfo = "text",
text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
type = "scatter"
) %>%
layout(
title = paste0(
"<b>Coverage PCA plot (<span style='color: blue;'>coverage > ",
coverageThreshold, "</span>)<b>"
),
xaxis = list(title = "<b><i>PC1<i><b>"),
yaxis = list(title = "<b><i>PC2<i><b>")
)
## Beta
pca.beta <- prcomp(t(beta))
p3.beta <- plot_ly(
x = pca.beta$x[, 1], y = pca.beta$x[, 2],
data = anno, color = ~Group, mode = "markers",
symbol = ~ as.numeric(as.factor(Group)),
marker = list(size = 10), hoverinfo = "text",
text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
type = "scatter"
) %>%
layout(
title = paste0(
"<b>Beta PCA plot (<span style='color: blue;'>coverage > ",
coverageThreshold, "</span>)<b>"
),
xaxis = list(title = "<b><i>PC1<i><b>"),
yaxis = list(title = "<b><i>PC2<i><b>")
)
# MDS plot
## Coverage
mds.cov <- cmdscale(dist(t(cov)))
p4.cov <- plot_ly(
x = mds.cov[, 1], y = mds.cov[, 2],
data = anno, color = ~Group, mode = "markers",
symbol = ~ as.numeric(as.factor(Group)),
marker = list(size = 10), hoverinfo = "text",
text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
type = "scatter"
) %>%
layout(
title = paste0(
"<b>Coverage MDS plot (<span style='color: blue;'>coverage > ",
coverageThreshold, "</span>)<b>"
),
xaxis = list(title = "<b><i>MDS1<i><b>"),
yaxis = list(title = "<b><i>MDS2<i><b>")
)
## Beta
mds.beta <- cmdscale(dist(t(beta)))
p4.beta <- plot_ly(
x = mds.beta[, 1], y = mds.beta[, 2],
data = anno, color = ~Group, mode = "markers",
symbol = ~ as.numeric(as.factor(Group)),
marker = list(size = 10), hoverinfo = "text",
text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
type = "scatter"
) %>%
layout(
title = paste0(
"<b>Beta MDS plot (<span style='color: blue;'>coverage > ",
coverageThreshold, "</span>)<b>"
),
xaxis = list(title = "<b><i>MDS1<i><b>"),
yaxis = list(title = "<b><i>MDS2<i><b>")
)
# SVD
## Coverage
svd.cov <- svd(x = t(cov))
p5.cov <- plot_ly(
x = svd.cov$u[, 1], y = svd.cov$u[, 2],
data = anno, color = ~Group, mode = "markers",
symbol = ~ as.numeric(as.factor(Group)),
marker = list(size = 10), hoverinfo = "text",
text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
type = "scatter"
) %>%
layout(
title = paste0(
"<b>Coverage SVD plot (<span style='color: blue;'>coverage > ",
coverageThreshold, "</span>)<b>"
),
xaxis = list(title = "<b><i>SVD1<i><b>"),
yaxis = list(title = "<b><i>SVD2<i><b>")
)
## Beta
svd.beta <- svd(x = t(beta))
p5.beta <- plot_ly(
x = svd.beta$u[, 1], y = svd.beta$u[, 2],
data = anno, color = ~Group, mode = "markers",
symbol = ~ as.numeric(as.factor(Group)),
marker = list(size = 10), hoverinfo = "text",
text = ~ paste("<b>Group:</b> ", Group, "<br><b>SampleID:</b> ", names),
type = "scatter"
) %>%
layout(
title = paste0(
"<b>Beta SVD plot (<span style='color: blue;'>coverage > ",
coverageThreshold, "</span>)<b>"
),
xaxis = list(title = "<b><i>SVD1<i><b>"),
yaxis = list(title = "<b><i>SVD2<i><b>")
)
# Clustering
## Coverage
p6.cov <- ggdendrogram(hclust(dist(t(cov)))) +
ggtitle(paste0("Coverage dendrogram (coverage > ", coverageThreshold, ")"))
## Beta
p6.beta <- ggdendrogram(hclust(dist(t(beta)))) +
ggtitle(paste0("Beta dendrogram (coverage > ", coverageThreshold, ")"))
list <- list(
cov, reshape.cov.melt, p2, beta, reshape.beta.melt, p1,
p3.cov, p4.cov, p5.cov, p6.cov, p3.beta, p4.beta, p5.beta, p6.beta
)
names(list) <- c(
"Coverage", "reshapedCoverage", "DHPlotCov", "BetaValues", "reshapedBetaValues", "DHplotBeta",
"PCAcov", "MDScov", "SVDcov", "HCcov", "PCAbeta", "MDSbeta", "SVDbeta", "HCbeta"
)
return(list)
}
Clustering of samples
Negative control
Simulated data
SessionInfo
## ─ Session info ──────────────────────────────────────────────────────────
## setting value
## version R version 3.5.1 (2018-07-02)
## os Ubuntu 16.04.5 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/Zurich
## date 2019-01-08
##
## ─ Packages ──────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 3.5.1)
## ade4 1.7-13 2018-08-31 [1] CRAN (R 3.5.1)
## ape 5.2 2018-09-24 [1] CRAN (R 3.5.1)
## assertthat 0.2.0 2017-04-11 [1] CRAN (R 3.5.1)
## autoplotly * 0.1.2 2018-04-21 [1] CRAN (R 3.5.1)
## backports 1.1.3 2018-12-14 [1] CRAN (R 3.5.1)
## bindr 0.1.1 2018-03-13 [1] CRAN (R 3.5.1)
## bindrcpp * 0.2.2 2018-03-29 [1] CRAN (R 3.5.1)
## Cairo 1.5-9 2015-09-26 [1] CRAN (R 3.5.1)
## callr 3.1.1 2018-12-21 [1] CRAN (R 3.5.1)
## cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.1)
## cluster 2.0.7-1 2018-04-09 [1] CRAN (R 3.5.1)
## ClusterR * 1.1.7 2018-12-10 [1] CRAN (R 3.5.1)
## colorspace 1.3-2 2016-12-14 [1] CRAN (R 3.5.1)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1)
## crosstalk 1.0.0 2016-12-21 [1] CRAN (R 3.5.1)
## data.table * 1.11.8 2018-09-30 [1] CRAN (R 3.5.1)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.1)
## devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.1)
## digest 0.6.18 2018-10-10 [1] CRAN (R 3.5.1)
## dplyr 0.7.8 2018-11-10 [1] CRAN (R 3.5.1)
## edgeR * 3.22.5 2018-10-02 [1] Bioconductor
## evaluate 0.12 2018-10-09 [1] CRAN (R 3.5.1)
## FD 1.0-12 2014-08-19 [1] CRAN (R 3.5.1)
## fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.1)
## geometry 0.3-6 2015-09-09 [1] CRAN (R 3.5.1)
## ggdendro * 0.1-20 2016-04-27 [1] CRAN (R 3.5.1)
## ggfortify * 0.4.5 2018-05-26 [1] CRAN (R 3.5.1)
## ggplot2 * 3.1.0 2018-10-25 [1] CRAN (R 3.5.1)
## glue 1.3.0 2018-07-17 [1] CRAN (R 3.5.1)
## gmp 0.5-13.2 2018-07-14 [1] CRAN (R 3.5.1)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 3.5.1)
## gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.1)
## gtools * 3.8.1 2018-06-26 [1] CRAN (R 3.5.1)
## hms 0.4.2 2018-03-10 [1] CRAN (R 3.5.1)
## htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.1)
## htmlwidgets 1.3 2018-09-30 [1] CRAN (R 3.5.1)
## httpuv 1.4.5.1 2018-12-18 [1] CRAN (R 3.5.1)
## httr 1.4.0 2018-12-11 [1] CRAN (R 3.5.1)
## jpeg 0.1-8 2014-01-23 [1] CRAN (R 3.5.1)
## jsonlite 1.6 2018-12-07 [1] CRAN (R 3.5.1)
## knitr 1.21 2018-12-10 [1] CRAN (R 3.5.1)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.5.1)
## later 0.7.5 2018-09-18 [1] CRAN (R 3.5.1)
## lattice 0.20-38 2018-11-04 [1] CRAN (R 3.5.1)
## lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.1)
## limma * 3.36.5 2018-09-20 [1] Bioconductor
## locfit 1.5-9.1 2013-04-20 [1] CRAN (R 3.5.1)
## magic 1.5-9 2018-09-17 [1] CRAN (R 3.5.1)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1)
## MASS 7.3-51.1 2018-11-01 [1] CRAN (R 3.5.1)
## Matrix 1.2-15 2018-11-01 [1] CRAN (R 3.5.1)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.1)
## mgcv 1.8-26 2018-11-21 [1] CRAN (R 3.5.1)
## mime 0.6 2018-10-05 [1] CRAN (R 3.5.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1)
## nlme 3.1-137 2018-04-07 [1] CRAN (R 3.5.1)
## OpenImageR 1.1.3 2018-12-09 [1] CRAN (R 3.5.1)
## permute 0.9-4 2016-09-09 [1] CRAN (R 3.5.1)
## pillar 1.3.1 2018-12-15 [1] CRAN (R 3.5.1)
## pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.1)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.1)
## plotly * 4.8.0 2018-07-20 [1] CRAN (R 3.5.1)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1)
## png 0.1-7 2013-12-03 [1] CRAN (R 3.5.1)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.1)
## processx 3.2.1 2018-12-05 [1] CRAN (R 3.5.1)
## promises 1.0.1 2018-04-13 [1] CRAN (R 3.5.1)
## ps 1.3.0 2018-12-21 [1] CRAN (R 3.5.1)
## purrr 0.2.5 2018-05-29 [1] CRAN (R 3.5.1)
## R6 2.3.0 2018-10-04 [1] CRAN (R 3.5.1)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.5.1)
## Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.1)
## readr 1.3.1 2018-12-21 [1] CRAN (R 3.5.1)
## remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.1)
## reshape2 * 1.4.3 2017-12-11 [1] CRAN (R 3.5.1)
## rlang 0.3.0.1 2018-10-25 [1] CRAN (R 3.5.1)
## rmarkdown 1.11 2018-12-08 [1] CRAN (R 3.5.1)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.1)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.1)
## shiny 1.2.0 2018-11-02 [1] CRAN (R 3.5.1)
## stringi 1.2.4 2018-07-20 [1] CRAN (R 3.5.1)
## stringr 1.3.1 2018-05-10 [1] CRAN (R 3.5.1)
## testthat 2.0.1 2018-10-13 [1] CRAN (R 3.5.1)
## tibble 1.4.2 2018-01-22 [1] CRAN (R 3.5.1)
## tidyr 0.8.2 2018-10-28 [1] CRAN (R 3.5.1)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
## tiff 0.1-5 2013-09-04 [1] CRAN (R 3.5.1)
## usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.1)
## vegan 2.5-3 2018-10-25 [1] CRAN (R 3.5.1)
## viridis * 0.5.1 2018-03-29 [1] CRAN (R 3.5.1)
## viridisLite * 0.3.0 2018-02-01 [1] CRAN (R 3.5.1)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1)
## xfun 0.4 2018-10-23 [1] CRAN (R 3.5.1)
## xtable 1.8-3 2018-08-29 [1] CRAN (R 3.5.1)
## yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.1)
##
## [1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.5
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library